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Thomas W. Strganac 
ABSTRACT 

A method for predicting unsteady, subsonic aeroelastic responses 
has been developed. The technique accounts for aerodynamic nonlinear- 
ities associated with angles of attack, vortex-dominated flow, static 
deformations, and unsteady behavior. The angle of attack is limited 
only by the occurrence of stall or vortex bursting near the wing. The 
fluid and the wing together are treated as a single dynamical system, 
and the equations of motion for the structure and flowfield are inte- 
grated simultaneously and interactively in the time domain. The 
method employs an iterative scheme based on a predictor-corrector 
technique. The aerodynamic loads are computed by the general unsteady 
vortex-lattice method and are determined simultaneously with the mo- 
tion of the wing. Because the unsteady vortex-lattice method predicts 
the wake as part of the solution, the history of the motion is taken 
into account; hysteresis is predicted. Two models are used to demon- 
strate the technique: a rigid wing on an elastic support experiencing 
plunge and pitch about the elastic axis, and an elastic wing rigidly 
supported at the root chord experiencing spanwise bending and twist- 
ing. The method can be readily extended to account for structural 
nonlinearities and/or substitute aerodynamic load models. The time 
domain solution coupled with the unsteady vortex-lattice method 
provides the capability of graphically depicting wing and wake motion. 
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CHAPTER I 


INTRODUCTION 

1.1 General 

Aeroelastic instabilities, which may be catastrophic, occur as a 
result of coupling between aerodynamic forces, structural forces, and 
mass. Therefore, the aeroelastic design of aircraft must address the 
structure and the aerodynamic forces within the flight environment. 
This concern includes the design of wings or aircraft for both flight 
and wind-tunnel testing. 

A method to predict the aeroelastic behavior of a wing is de- 
scribed. The objective is to develop an improvement in aeroelastic 
analyses by accounting for aerodynamic nonlinearities associated with 
angles of attack, static deformations, vorticity-dominated flow, and 
unsteady behavior. This method generates a realistic simulation of 
aeroelastic response or instability by predicting both the motion of 
the wing and the motion of the fluid simultaneously. In other words, 
the wing and fluid are treated as a single dynamical system, and the 
equations of motion for the structure and flowfield are integrated 
simultaneously and interactively in the time domain. 

The equations of motion are developed for classical two-degree- 
of-freedom wing motion; that is, a rigid wing which is allowed to 
plunge and pitch about the elastic axis. In addition, the equations 
of motion are developed for an elastic wing with bending and twisting. 
The unsteady vortex lattice method (UVLM) is used to predict the 
aerodynamic loads. The technique accounts for the aerodynamic 
nonlinearities. 
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A complication does exist with a time domain approach: the aero- 

dynamic loads cannot be predicted unless the motion of the wing is 
known, and the motion of the wing cannot be predicted unless the aero- 
dynamic loads are known. An iterative scheme that accounts for the 
interaction between the aerodynamic loads and wing motion was devel- 
oped based upon a predictor-corrector approach. An advantage of the 
method is that it provides an excellent opportunity to animate the 
motions in the flowfield and wing since the governing equations are 
solved in the time domain. 

The technique described herein provides a somewhat different ap- 
proach to model flutter for a wing of arbitrary planform. The equa- 
tions are developed about arbitrary static angles of attack, providing 
a new capability to study the associated nonlinear effects. Through 
the use of the UVLM one is able to treat low-aspect-ratio wings and 
account for the wake rollup at the wing tips. The UVLM models the 
wake where the history of the motion resides. 

1.2 The Aeroelastic Phenomenon 

Many different types of aeroelastic instabilities exist. Static 
instabilities include aeroelastic divergence in which the elastic re- 
storing forces of the wing are exceeded by the aerodynamic forces. 
Dynamic instabilities such as flutter are caused by an exchange of 
energy between the aircraft and surrounding air. General aeroelastic 
phenomena are described by Fung (1955), Bisplinghof f , Ashley and 
Halfman (1955), and Dowell (1980) among others. A historical per- 
spective of aeroelastic research is given by Collar (1978). 
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A complete model of the aerodynamic forces, elastic forces, and 
inertial forces is required to describe the unsteady aeroelastic na- 
ture of an aircraft. As Figure 1.1 depicts, this interaction among 
these three types of force provides the conditions necessary for 
potential static and dynamic instabilities. The nature of the aero- 
elastic problem can best be illustrated by this three-ring diagram as 
suggested by Yates (1971). Variations to this diagram are presented 
by Bisplinghoff et al. (1955), Collar (1978), and Dowell (1980). 

Aerodynamic, elastic, and inertial forces are shown. The common 
ground between the inertial forces and elastic forces represents the 
f ree-vibration problem, the common ground between the aerodynamic 
forces and elastic forces represents static aeroelastic phenomena 
which include instabilities such as divergence, and the common ground 
between the aerodynamic forces and inertial forces represents dynamic 
stability investigations of "rigid" aircraft. The interaction of all 
three forces is present in the aeroelastic phenomenon known as 
flutter. 

Flutter has been defined in many ways. Among the first to exam- 
ine the nature of flutter, Theodorsen (1940) described flutter and 
developed a general theory which is still used today. Pines (1958) 
defines flutter as "a self-excited oscillatory aerodynamic instabil- 
ity" and states that flutter may be considered "an oscillatory in- 
stability where one degree of freedom is driven at resonance by a 
second degree of freedom, both at the same frequency". Hancock, 
Wright, and Simpson (1985) describe flutter as "a complex phenomenon 
where, in the classical sense of the term, two or more structural 
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normal modes are coupled and excited through time dependent aero- 
dynamic loads". Dowell (1980) also describes the flutter phenomenon 
and suggests that many types of flutter exist; these include: 
"'Coalescence' or 'merging frequency' flutter", " 'Single-degree-of- 
freedom' flutter", and "'Divergence' or 'zero frequency' flutter". 

An illustration of the flutter phenomenon is presented in Fig- 
ure 1.2. A wing is initially moving through a fluid at rest. The 
wing is subjected to a disturbance. In the left-hand side of the 
figure the dynamic pressure is below the critical condition for flut- 
ter. The bending and torsion of the wing occur at separate, unique 
frequencies which are close to the natural frequencies for the free 
vibration modes. The motion decays following the initial disturbance 
as there is not enough energy being extracted from the freestream to 
maintain the motion. In fact, the freestream is extracting energy 
from the wing and damping the motion. There is aerodynamic damping 
here. The system is being damped by aerodynamic effects only. In the 
right-hand side of the figure the dynamic pressure is greater than the 
critical pressure. Not only is the motion seen to grow in this case 
but bending and torsion occur at a common frequency - a characteristic 
associated with flutter. In the lower portion of the figure the shift 
in system frequencies due to the dynamic pressure is illustrated. The 
frequencies are the natural ones when the dynamic pressure is zero. A 
coalescence occurs at a critical dynamic pressure, the system is now 


in flutter 



PLUNGE PITCH 
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1.3 Literature Review 

The purpose of this research is to provide an improvement in 
aeroelastic analysis capabilities. Through the application of the 
UVLM to the aeroelastic formulation new capabilities are introduced. 

As mentioned previously these include the modeling of aerodynamic 
nonlinearities and new solution strategies. These enhancements, in 
turn, open the door for further improvements primarily in structural 
nonlinearities. However, the emphasis in this research is placed on 
the aerodynamic model and the solution strategy; the literature review 
will focus on these areas. 

1.3.1 Aeroelastic models 

Aeroelastic research is very broad indeed. Experimental and 
theoretical research has been performed in structures, fluids, and 
dynamics as well as solution schemes associated with aeroelasticity. 

Desmarais and Bennett (1978) described an approach for flutter 
analysis which is representative of current numerical strategies. 

In their method the generalized mass and mode shapes are developed 
through an external structural analysis package. These modes are used 
to describe the shape of the wing which in turn is used to compute an 
aerodynamic influence matrix. A subsonic kernal function approach is 
used to model the aerodynamics. An eigenvalue formulation is used 
with the governing equations. Hence, the solution is performed in the 
frequency domain, the governing equations are cast in a linear form, 
and a simple harmonic form of the solution is assumed. As we shall 
see, the method developed in this research is not limited by these 
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assumptions. The roots of this eigensolution, which are dependent 
upon the velocity, density, and Mach number, are used to determine the 
critical conditions necessary for flutter. 

The method commonly used to determine the flutter conditions of a 
wing using a frequency domain solution is illustrated in Figure 1.3. 
The method, called the V-g (for Velocity-incremental damping) method, 
tracks the roots of the modes used for the solution for increasing 
values of velocity. A fictitious damping (the root) is determined at 
each velocity which is required to satisfy the eigenvalue problem. In 
the upper part of the figure we plot the roots for a 3 mode system. 

The upper plane represents an unstable region; hence, when a root 
enters this region the system is unstable. At only this "flutter 
crossing" is the solution meaningful, yet no physical interpretation 
of the motion is available. In the lower part of the figure the 
frequencies for these modes are computed as the velocity increases. 

Hassig (1971) and Fung (1955) described some alternate solution 
schemes for the frequency domain solutions. Fung also presented many 
examples, some of which will be referred to in this dissertation. 
Goland and Luke (1949) outlined a method to describe the aeroelastic 
character of an elastic wing. 

Garrick (1969) presented several classical papers on the flutter 
phenomenon, flutter determination, aerodynamic theories, and other 
associated subjects. Loring (1941) introduced generalized coordinates 
through the use of the natural modes to solve the flutter problem. 

His work developed a systemized approach for the structural model. 



FREQUENCY INCREMENTAL DAMPING 
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Figure 1.3. Determination of Flutter for Frequency Domain 
Formulations, V-g and Frequency Diagrams. 
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Dowell (1980) described both frequency and time domain solutions 
of the aeroelastic equations* Frequency domain solutions are the most 
widely used since these require less resources and a more convenient 
solution scheme. However, a frequency domain solution does not pro- 
vide a physical description of the motion as do time domain solutions. 
Time domain solutions, which do not assume a form of the solution 
(i.e. harmonic motion) as do the frequency domain schemes, provide the 
more general and physical approach. 

A more complete formulation of the aeroelastic equations was 
described by Bisplinghoff (1975) and Fung (1955). Descriptions of 
aeroelastic phenomena were presented by Theodorsen (1940), Biot and 
Arnold (1948), Pines (1958), and Hancock et al. (1985). These later 
two papers presented physical models of flutter. Pines, in par- 
ticular, examined the governing equations in a general sense and 
established relationships for which flutter can exist. 

An aeroelastic model, which is similar in nature to the one pre- 
sented in this dissertation, is described by Devers (1972). Devers 
developed a time domain solution of the governing equations, he too 
used a vortex-lattice method as an aerodynamic model. However, sim- 
plifications in his model prevented the unsteady nature of trailing 
edge and wing-tip vortex effects to be considered. In addition, his 
integration scheme assumed a form of the solution and, as a conse- 
quence, the simulation reflected this assumption. 
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1.3.2 Aerodynamic models 

Currently, the major effort of aeroelastlc research is to de- 
scribe the unsteady aerodynamic loads, and, in particular, the flow in 
the transonic flow regime as most recently described by Batina et al. 
(1987). Edwards (1986) presented an overview of aerodynamics associ- 
ated with transonic aeroelasticity. 

Bisplinghoff (1975) described many different aerodynamic opera- 
tors which are available for all flow regimes. Yates (1985) described 
the various approaches available for aerodynamic models, which include 
doublet-lattice, kernal function, velocity potential and acceleration 
(pressure) potential formulations. 

The research will concentrate in the subsonic regime. Theodorsen 
(1940) introduced a two-dimensional functional form of the subsonic 
aerodynamic operator which is still popular today. Desmarais and 
Bennett (1978) implemented a subsonic kernal function indicative of 
more current strategies. However, most methods are quasi-steady, lin- 
ear (i.e. small angle) theories. A new approach to the subsonic aero- 
dynamic model for aeroelastlc analysis is introduced. The unsteady 
vortex-lattice method (UVLM) provides the opportunity to address aero- 
dynamic nonlinear effects associated with unsteady flow aspect ratios, 
static deformations, and the angles of attack. Hence, nonlinear 
aerodynamic effects may be studied as they apply to the aeroelastic 
phenomenon. 

Experimental investigations have demonstrated nonlinear effects. 
Yates and Bennett (1971), Farmer (1982), and Yates, Wynne, and Farmer 
(1982) addressed the effect angle of attack has on flutter boundaries. 
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Hsu (1957) addressed the effect aspect ratio has on flutter 
boundaries * 

The development of the UVLM will be briefly reviewed. The reader 
is referred to the work of Kandil (1974), Maddox (1973), and most re- 
cently, Konst ad inopoulos (1981) for more complete descriptions. 

The unsteady vortex-lattice method, as we use it in this research, 
originates from the investigation of Belotserkovskiy (1966). This 
original work could treat arbitrary wing planforms and deformations 
but could not model the geometry of the wake which thereby limited it 
to small angles of attack. Belotserkovskiy (1968) then developed a 
nonlinear method to calculate the aerodynamic loads on wings with 
wing-tip separation in steady flows. Later, Belotserkovskiy and Nisht 
(1974) presented a method for the treating of rather general planforms 
in nonlinear unsteady flow. As part of this research, they determined 
the shape of the wake convecting from the wing tips and the trailing 
edge. 

Maddox (1973) and Mook and Maddox (1974) used this vortex-lattice 
method and considered leading edge separation. However, this method 
did not account for force-free wing-tip and trailing-edge vortex 
sheets. Kandil (1974) and Kandil, Mook, and Nayfeh (1974) refined 
this approach and their results were in very good agreement with the 
experiments. The remarkable accuracy of vortex-lattice methods was 
discussed by James (1971). 

Atta (1976) and Atta, Kandil, Mook, and Nayfeh (1976, 1977) ex- 
tended the method to treat unsteady flows past rectangular wings with 
sharp-edge separation. The problem was posed in an inertial frame 
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which proved to be awkward. Thrasher, Mook, Kandil, and Nayfeh (1977) 
and Thrasher (1979) posed the problem in terms of a body-fixed frame. 
The approach treated rectangular wings executing arbitrary maneuvers 
as long as separation occurred along the wing’s sharp edges and vortex 
bursting did not occur at or near the wing's surface. 

Thrasher (1979) coupled the aerodynamic method with a predictor- 
corrector method to predict the flowfield, loads, and motion of a 
hinged rectangular wing due to a prescribed motion of a flap. This 
development is of particular interest since it had application to 
aeroelastic type problems. Kandil, Atta, and Nayfeh (1977) and Atta 
(1978) refined the approach to treat delta wings. Further, they 
developed a higher-order convection process for the wake, but this 
process required iteration. 

Nayfeh et al. (1979) modified the method to treat small, harmonic 
oscillations around an arbitrary angle of attack. The general 
unsteady method was described by Nayfeh, Mook, and Yen (1979), 
Konstadinopoulos (1981) and Konstadinopoulos, Mook, and Nayfeh (1981). 
Most recently the general method was described by Konstadinopoulos, 
Thrasher, Mook, Nayfeh, and Watson (1985). 

Recently, Konstadinopoulos (1984) and Konstadinopoulos, Mook, and 
Nayfeh (1985) developed a numerical method for simulating subsonic 
wing rock. More recently, Elzebda (1986) described a model for sim- 
ulating two-degree-of-f reedom wing rock. Elzebda, Mook, and Nayfeh 
(1985) described the ability of the general method to model the aero- 
dynamics for close-coupled lifting surfaces. 
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1.3.3 Structural models 

In this research the wing structure is represented by two differ- 
ent models for demonstration purposes - a rigid wing which can pitch 
and plunge about the elastic axis and an elastic wing. Currently, for 
the model of the elastic wing, the elastic axis of the wing is repre- 
sented by a beam which is allowed to bend and twist. No structural 
nonlinearities are taken into account. The two-degree-of-f reedom 
model is addressed in the texts on aeroelasticity we have already 
cited. Several texts (for example, see Rivello (1969)) describe the 
equations of a beam which we have implemented. 

The numerical model can be readily extended to include nonlinear 
effects of beam flexure and torsion. Several references consider 
nonlinear models, these include the work of Crespo da Silva and Hodges 
(1986), Hodges and Dowell (1974), and Kaza and Kvaternik (1977). 

These references address the more difficult kinematics and large 
deformations associated with rotorcraft applications. 

More complicated models of the wing are typically handled by fi- 
nite element methods. For example, the method described by Desmarais 
and Bennett (1978) implements the structural modes generated by any 
experimental or analytical method. Hence, the description of struc- 
tures is dependent on the capability of these alternate sources. Yet, 
an advantage in the aeroelastic design of structures is gained through 
the inclusion of more elaborate structural models. Ashley et al. 
(1980), Bendiksen and White (1986), and Weisshaar (1987) describe the 
advantages of aeroelastic tailoring of structures. 
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1.4 The Present Method 

A method for predicting unsteady, subsonic aeroelastic responses 
is described. Previously, we described the application of the UVLM to 
the two-degree-of-f reedom problem (Strganac and Mook (1986)). The 
introduction of the UVLM to the prediction of flutter provides the 
opportunity to model aerodynamic nonlinearities. Later, we extended 
this effort to include the elastic wing model (Strganac and Mook 
(1987)). Most recently, we demonstrated the simulation of flutter by 
animating the wing and flowfield as predicted by our model (Strganac, 
Mook, and Mitchum (1987)), These efforts are further described here. 

In Chapter II, the equations of motion for a rigid wing on an 
elastic foundation are developed. Structural damping and static 
deflections are included. The nondimensionalization concept is also 
introduced. 

In Chapter III, the equations of motion for an elastic wing on a 
fixed support are developed. The equations account for static defor- 
mations. Mass and stiffness matrices for the wing are developed. 
Coupling exists between bending and torsion. Structural nonlineari- 
ties are not addressed but the general formulation allows for this 
feature. 

In Chapter IV, the unsteady vortex-lattice method is described. 
The computation of the wake geometry is illustrated. Comparisons of 
the computed pressure distribution and the sensitivity to aspect ratio 
are shown. 

In Chapter V, the integration schemes which account for the 
structural dynamic and aerodynamic interaction are developed. The 
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methods are based upon a predictor-corrector technique, with appro- 
priate modifications to account for aerodynamic/dynamic coupling. The 
convection schemes which are used to accurately generate the wake in 
the unsteady model are discussed. The rigid wing and elastic wing 
formulations are individually described as each present unique 
requirements. 

In Chapter VI, several examples of computed results are pre- 
sented. Wing response, aerodynamic loads, and pressure distributions 
are shown. In addition, phase planes and aerodynamic load histories 
are shown which further describe the physical nature of aeroelastic 
phenomena. Finally, the ability to take advantage of time domain 
solutions is graphically demonstrated. The wing and wake shape and 
vorticity distribution are shown at several time steps. These 
sequences are indicative of the animation capability which exists 
with the present formulation. 



CHAPTER II 


EQUATIONS OF MOTION FOR THE RIGID WING 
2.1 Description of the System 

In this chapter the equations of motion are developed for a rigid 
wing on an elastic support limited to plunge and pitch about the 
elastic axis. This simplified structural model, illustrated in Fig- 
ure 2.1, serves as a demonstration of the technique and is described 
in many references (see Fung (1955)) as a common example. However, 
the equations are developed about an arbitrary static angle of attack, 
<J>, which is included in the aerodynamic model. 

The wing is attached to the support by linear springs to model 
stiffness characteristics and by linear dashpots to model damping 
characteristics. The point of attachment is referred to as the elas- 
tic axis. The center of gravity is offset from the elastic axis by 
the parameter r. This offset affects sensitivity to flutter; this 
effect has been described by Pines (1958) and Bisplinghoff (1975), as 
well as others. The aerodynamic model accounts for spanwise effects; 
hence, the aspect ratio of the wing is considered. Physical proper- 
ties of the wing, stiffness (and structural damping) properties of 
the support, and aerodynamic properties act on the wing section as 
illustrated. 
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Figure 2.1. The Rigid-Wing Model. 
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2.2 Development of the Equations 

The components of the acceleration of the center of gravity are 
described in the inertial frame. The x-component is expressed as 


-ra sin(a +<(>)- rot cos(a + 


( 2 . 1 ) 


and the y-component is expressed as 

A 

a y “ y + r ® cos(a +<)>)- ra sin(a + <(») (2.2) 

where a and y are the coordinates which describe the pitch and 
punge motion. The location of the center of gravity relative to the 
elastic axis is given by r, which is a positive quantity if the 
center of gravity is aft (i.e., towards the trailing edge) of the 
elastic axis. 

Summing forces in the X direction results in the following 
equation: 

F x + N sin(a + 4> ) 

“ m{-ar sin(a + <(>)- a 2 r cos(a + <(»)} (2.3) 

Summing forces in the Y direction results in the following equation: 

Fy - N cos(a + $) 

•• •• .2 

* m{y + ar cos(a +<)>)- a r sin(a + <j>)} 


( 2 . 4 ) 
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Summing moments about the wing’s center of gravity results in the 
following equation: 


M r + M g + Nx - rFy cos (a + <j>) 

+ r sin(a + <j>)F * I a 
x eg 


(2.5) 


where x is the distance between the reference position of the aero- 
dynamic loads and the mass center of the wing and is a positive 
quantity if the mass center is aft of the reference position. The 

reactions of the springs and dashpots are denoted by F F , and 

x y 

M g . Equations (2.3) and (2.4) are substituted into Equation (2.5), 
and the result is the following differential equation: 


I a = Nx + M + M - rN - mry cos(a + 6) - ma 
eg r s 


( 2 . 6 ) 


Equations (2.4) and (2.6) represent the differential equations 
which describe the motion of the system. Linear elastic spring and 
damping forces and moments are defined as follows: 

. 

F = -c _y - k y 

y T y 


These expressions are substituted into Equations (2.4) and (2.6). 
Structural nonlinearities could be introduced in the spring and 


damping model 
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The force and moment equations are now rewritten as 


«• •• 


my + mra cos(a + $) - mo r sin(a + <f>) 


-N cos (a + <)>)- k^y - C^y 


(2.7) 


I a + r ma + rmy cos (a + <j») 
eg 


• 2 •• 
Nx~ca-ka + M -rN-r ma 
a a r 


( 2 . 8 ) 


The form of these equations can be simplified by introducing the 
following: 


I ■ I + r m 
e eg 


and 


M - + (x - r)N 


Equations (2.7) and (2.8) are now rewritten as 


I a + mry cos (a + 4) + k a + c a = M 
e a a 


(2.9) 


and 


my + mra cos(a + <)>)- mra sin(a + 4>) 


+ ky + cy“-N cos (a + <p) 

y y 


( 2 . 10 ) 


The linear form of Equations (2.9) and (2.10) agrees with the 
equations derived by Fung (1955). 
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The numerical integration scheme which is described in Chapter V 
requires first-order differential equations. To this end. Equa- 
tions (2.9) and (2.10) are solved simultaneously to provide equations 
in terms of the second-order time derivatives. These equations are: 


2 2 

, I - rar cos (a + <(>) 
L e 


{ 


(k a + c a)r cos(a + <fr) 
a a 


+ I o^r sin(a + <|>) — (k y + c y) 

e m y J y J 


- Mr cos(a + ) 



cos (a + $) 
m 


} 


( 2 . 11 ) 


2 2 
I - mr cos (a + <J> ) 


L e 


{-c a - k a + M 
1 a a 


2.2 

-ram sin(a + <f>) cos(a + 4*) 


A 

+ r cos (a + 4)(k^y + c^y) + Nr cos (a + 4)} 


( 2 . 12 ) 


where the aerodynamic normal force, N, and pitch moment, M, are 
defined as 


N 


1 2 
T P U « AC „ 

2. 00 n 


(2.13) 


M = i pU^Ani C 
2 « cm 


(2.14) 


M is referenced to the elastic axis. The quantity n^ c is the full 
chord of the wing. 
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2.3 Nondimensionallzed Form of the Governing Equations 

The governing equations of the UVLM have been written in a non- 
dimensional form through the introduction of the characteristic 
length, fc c , characteristic velocity, and the resulting charac- 

teristic time, (x * £ C /Uoo)* Equation (2.11) is also rendered non- 

2 

dimensional by dividing both sides of the equation by U /A . This 

oo C 

equation is now expressed as 


Y 


1 

K - r^ cos^(a + $) 

A 

+ K{a r sin(a + <J>) 


[<V + 

- K Y - 

y 


C a - 
a 


C Y - 

y 


nCC )r cos (a + $) 
m 

CC cos(a + $) } ] 
n 


( 2 • 15) 


2 

where K ■ I /m£ is introduced as a nondimensional inertia and 
e c 

C * A/2m is introduced as a nondimensional aerodynamic constant# 

Hence, we also utilize the mass, m, to characterize the physical 

properties in nondimensional form# Other terms in Equation (2.15) 

have been rearranged and represent the necessary collection of terms 

for the nondimensional form. All quantities are nondimensional. 

Similarly, Equation (2.12) is nondimensionalized by dividing both 

2 2 

sides of the equation by This equation is now expressed as 

a — L_ {(C y + K Y)r cos(a + +) 

K - r cos (a + <|>) ^ Y 

A A ^ 

- r a sin(a + <(> ) cos(a + <J>) - (C a + K a) 


+ C(n + r cos^(a + <J>)C n )} 


(2.16) 
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One should note that the velocity now appears only in the spring and 
damping terms. K 0 and Ky contain the square of the speed in the 
denominator. C a and contain the speed in the denominator. 

Hence, the effective stiffness and damping of the system decreases as 
the freestream speed increases. 



CHAPTER III 


EQUATIONS OF MOTION FOR THE ELASTIC WING 

3.1 Description of the System 

In this chapter the equations of motion for an elastic wing are 
developed by using an energy approach. This wing is illustrated in 
Figure 3.1. Small displacements about an equilibrium position are 
assumed. The angle of attack is introduced into the equations. The 
wing is allowed to bend and twist about an elastic axis not necessar- 
ily coincident with the inertial axis. Thus, there is inertial as 
well as aerodynamic coupling between wing bending and torsion. 

The wing is represented as a cantilevered beam, which can both 
bend and twist. The cross section of the beam is assumed to be rigid. 
The wing root is rigidly fixed. The physical properties may vary 
along the span of the wing. Spanwise variation in the aerodynamic 
loading is considered. The solution of the governing equations con- 
sists of static and dynamic contributions. The statically deformed 
shape due to steady aerodynamic loads and the distributed weight is 
determined. Franz, Krenz, and Kotschote (1984) discussed the impor- 
tance of deformations due to static loads on aerospace vehicles, 
particularly wind-tunnel models. 

3.2 Development of the Equations 

We ignore the out-of-plane (that is, the spanwise) motion. The 
components of the velocity are expressed as 

v x - -x a (y) sin(a(y,t)) a(y,t) (3.1) 
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Figure 3.1. The Elastic-Wing Model 



The kinetic energy of the system is expressed as 


T * / i [m(y)({w + x (y)(cos a)a} 2 
0 

+ { x (y ) ( s in a)a} 2 ) + I (y)a 2 ] dy (3.3) 

eg 

The potential energy of the system is expressed as 

V - / {EI(y)w" 2 + GJ(y)a' 2 - m(y)g(w + x (y) sin o) } dy (3.4) 

0 a 

The distributed aerodynamic loads are provided by the unsteady 
vortex-lattice method (UVLM). The work done by the aerodynamic load- 
ing is expressed as 

L 

W * / (-Nw(y,t) + Mct(y ,t) } dy (3.5) 

0 

This aerodynamic spanwise loading is dependent upon the motion, 
and the motion, in turn, is dependent upon the aerodynamic loads. 

The normal force and pitching moment are functions which should be 
expressed as 

N = [w, w, w' , a, a, t; (j), AR, planform] 

M * [w, w, w' , a, u, t; 4>, AR, planform] 
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According to Hamilton's principle, 
t L 

fl (6T - 6V + 6W) dy dt = 0 (3.6) 

0 0 

Substituting T, V, and W into Equation (3,6), one obtains the 
following differential equations: 

*• *2 
mw + mx cos(a)a - mx sin(a)a 
a a 

+ (El w")" - mg - -N (3.7) 

•• •• 

I a + mx cos(a)w - (GJ a’) 1 “ mgx cos(a) * M (3.8) 

e a a 

The linear form of these equations is 

mw + mx a + (El w") tf - mg = -N (3.9) 

a 

I^a + mx^w - (GJ a 1 )’ - mgx^ * ^ (3.10) 

where all properties vary along the length. 

For the case where Equations (3.9) and (3.10) represent a canti- 
levered wing, the boundary conditions for these equations are 

w(0,t) * w T (0,t) * w ,? (L, t ) * (EI(L)w ff (L,t)) f * 0 
and (3.11) 

a(0, t) * a f (L, t) 38 0 
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3.3 Static and Dynamic Equations 

The solution of Equations (3.9) and (3.10) provides both the 
static and dynamic contributions. The developed equations will yield 
the static solution for cases of decaying motion. However, the solu- 
tion may require unnecessary computations for static solutions only; 
therefore, w(y,t) and a(y,t) are redefined as 

w(y,t) - w g (y) + w^(y,t) (3.12) 

and 

a(y,t) - a g (y) + a d (y,t) (3.13) 

where the subscripts s and d refer to the static and dynamic solu- 
tions. Now these equations become 

mw . + mx a . + (EIw")" + (EIw")" - mg - -N (3.14) 

dad d s 

and 

I a + mx w . - (GJa')' - (GJal)' - mgx * M (3.15) 

e d a d s d a 

From Equations (3.14) and (3.15), the equations for the static 
deformations and the dynamical equations for small motion about the 
static deformations are 

(Elw'g)" - mg - -N 

-(GJa;)* - mgx a - M g 

raw. + mx a, + (EIw'J)" ■ -N , 
a ot a a a 


(3.16) 

(3.17) 


(3.18) 
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Vd + "Vd " (GJa d r " M d (3 * 19) 

Meirovitch (1980) describes a method of solution to these equations by 
an expansion of the dependent variables* The variables w and a 
are represented by expansions in terms of the natural free vibration 
modes of the system* 


I 

w * 51 q. (3.20) 

i=l 1 


and 


a 


J 


I 

j-i+i 




(3.21) 


where the and qj are the generalized coefficients and the 

and are the modes chosen as the comparison functions, and where 

I * number of bending modes selected, 

J - I = number of torsion modes selected, and 

J ■ total number of modes selected to represent the solution. 


3.3.1 Equations governing the static contribution 

These expansions are substituted into the static Equations (3.16) 
and (3.17). The result becomes 

I 

l (El*;)" q. - mg - -N (3.22) 

i-i 1 1 s 


J 

l (GJ*’) f q. ■ M + mgx 
j-I+1 J J S 


(3.23) 
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Multiplying the first equation by * r for r=l,2,...,I and then 
integrating along the span, one obtains 


LI L 

/ £ dy q, - / mg* dy = -/ N $ dy 

0 i=l 1 r 1 0 r s r 


(3.24) 


Multiplying the second equation by $ g for S=I+1,...,J and 
then integrating along the span, one obtains 


L J 


*/ l (GJ*!)’* dy q 
0 j=I+l J J 


L L 

/ M * d y + / ®g x a * s d y 

0 S s 0 as 


(3.25) 


In Equations (3.24) and (3.25), the first term is integrated by parts 
and the boundary conditions, Equation (3.11), are imposed. With rear- 
rangement these equations are rewritten as 


I L L L 

I J EI*V*" dy q = / mg* dy - / N * dy 

_ i a ir 1 a r a s r 


(3.26) 


i-1 0 


O 


and 


J L 

I I GJ*!** dy q 

j-I+1 0 J 3 


L L 

f M * s dy + / mgx * dy 
0 s s Q as 


(3.27) 


The above expressions represent J equations, which in matrix 
form may be written as 


[K] {q} = [{A }. + {S}] 

s s 


(3.28) 
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The stiffness matrix [K] is partitioned as 


K 11 K 12 


L K 12 K 22 J 


where 


c n - / 


for i * 1,2, ... ,1 

j “ 1,2,...,! 


K 12 ' K 21 ’ 0 


t 22 - J dy 


for 1. 88 1^* 1 y • • • I «J 

j 88 1+1 , • • • , J 


{A g } is the static aerodynamic loading vector and is defined as 


r l 

/ - Vi 

0 s 1 


<V » < 


> 


f M * dy 
J sj 


VO 


for i=l ,2, . . . ,1 
J “H" 1 , • • • , J 


as 


{S} is the static loading (due to weight) vector and is defined 


C L 

/ mg* dy 

A ^ 




{S} - < 


> 


/ mgx„ * dy 


VO 


a j 


J 


for i=l ,2, . . . ,1 
j =1+1 , . . . , J 
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3.3.2 Equations governing the dynamic contribution 

We now concentrate on the dynamic equations. Substituting the 
series representation into the governing equations ((3.18) and (3.19)) 
for the dynamic solution yields the following equations: 


I 

l 


J 

1 


mq ♦ + l mx q * + 

i-1 j-I+1 3 3 


I 

l 

i-1 


q i (EI*p'’ = -N d 


I 

I 


mx q . + 
i-i a * 1 j.i + i 


j 

i 


Wj ' j«x+l 


M, 


(3.29) 


(3.30) 


Multiplying the first equation by $ r for r - 1,2,..., I and 
integrating along the span we obtain the following expression: 



I L L 

+ l q. ! dy = -/ N * dy (3.31) 

i-1 0 1 r 0 a r 


Multiplying the second equation by $ g for s - R+1,...J and 
integrating along the span we obtain the following expression: 


I 

l 

i-1 


mx ♦. # 
a i s 


dy + 


J 

I 

j+I+1 


L 

J I 
0 


e j s 


dy 


+ I 

j-I+1 


q j 


L 

J GJ*!*’ dy 



M, ♦ dy 
a s 


(3.32) 
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The third term in each of these expressions has been integrated 
by parts and the boundary conditions, Equation (3.11), are imposed. 

The result is a more simplified form of the equations. The above 
expressions represent a number of equations (J) equal to the total 
number of modes chosen to represent the system. These equations may 
be written in the following matrix form: 

[M]{qJ + [K]{qJ - {A^ (3.33) 


The mass matrix [M] is partitioned as 


[M] 


" M ll 

*12“ 

_ M 21 

1 

CM 

CM 

X 


where the terms in [M] are defined as follows: 


L 

M , , ■ / m#. ♦. dy 

11 0 1 J 


for i 

j 


1.2. . ...1 

1.2.. ..,! 


M 


12 


M, 


21 



mx ♦. ♦, 
a 1 j 


dy 


for i * 1,2, ... ,1 
j ™ 1+1 »••*»•! 


M, 


22 


L 

J I 
0 


*, *. 
e i j 


dy 


f or i * 1+1 , . . . , J 
j ■ I+1,...,J 


The stiffness matrix [K] is as previously defined 
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The right-hand side contains the aerodynamic contribution and is 
expressed as follows: 


I -Vi dy 

0 a 1 


<v 




J M # dy 
^0 d J ^ 


for i ■ 1,2, ... ,1 


and represent the sectional normal force and moment refer- 

enced at the elastic axis and computed from the dynamic solutions. 

[M] and [K] are symmetric. In addition, these equations may 
be inertially uncoupled if the sectional center of gravity and elastic 
axis are coincident such that x Q * 0. However, the equations are 
still coupled through the aerodynamic loads. 

The equations are to be cast in a manner similar to that of the 
two-degree-of-f reedom formulation. Therefore, n * q is introduced 
which allows the governing equations to be written in state space 
form. 


f{n} 1 

’ 0 — [M ] ~ 1 [K]' 

f{r>} 1 


"[M]" 1 {A}/ 

J t » 


\ 

+ 

r 

y q} d\ 

.[I] o 

Kj 

1 

l o J 


The above expressions provide a set of first-order differential 
equations, the solution of which describes the motion of the elastic 
wing. The integration scheme of these equations will be described 
subsequent ially. 



CHAPTER IV 


AERODYNAMIC LOADS, THE UNSTEADY VORTEX-LATTICE METHOD 
4.1 Overview 

The equations of motion and the associated integration scheme are 
formulated in a manner which permit the aerodynamic loads to be deter- 
mined concurrently with the motion. In this chapter, a technique is 
described which models the aerodynamic loads using the unsteady vortex 
lattice method (UVLM). This method has been completely described by 
Konstadinopoulos (1981); hence, the technique as applied in this 
research will be briefly reviewed. 

The UVLM is a numerical model of the three-dimensional flowfield, 
which can be used to treat arbitrary maneuvers of wings of arbitrary 
planforms, including highly swept delta wings which exhibit leading 
edge separation, multiple closely coupled lifting surfaces, and low- 
aspect-ratio planforms. It can also treat arbitrary angles of attack 
and camber as long as stafl or vortex bursting in the near wake does 
not occur. The method accounts for the nonlinear effects of the wakes 
adjoining the tips and trailing edge. 

The ability of the UVLM to predict accurate spanwise distribu- 
tion of loads has been demonstrated. In addition to the work of 
Konstadinopoulous, it is also worth mentioning other applications of 
the UVLM. The method was used to model aerodynamics of close-coupled 
lifting surfaces (Elzebda (1986)). Mook and Nayfeh (1985) demon- 
strated the method for high-angle-of-attack aerodynamics. Kobayakawa 
and Onuma (1985) applied the vortex-lattice method to model propeller 


36 



37 


aerodynamics. Thomas and Nerney (1976) implemented the vortex-lattice 
method (coupled with slender body theory) to predict the aerodynamics 
of wing-body combinations. Kandil, Mook, and Nayfeh (1976) applied 
the UVLM to aircraft interference problems. 

The approach assumes that the flow is incompressible and inviscid 
and does not separate on the wing, but that separation occurs along 
the sharp edges where the Kutta Condition is imposed in a steady flow. 
At each time step this vorticity, which forms the wake, is convected 
at the local particle velocity; thus, the position of and the distri- 
bution of vorticity in the wake are part of the solution. The wake 
contains the history of the flow. Consequently, the velocity induced 
by it at the present time reflects the previous motion. 

In the aerodynamic model the vorticity in the flow is restricted 
to a thin region around the lifting surface and its wake. The flow 
outside this region is considered irrotational. The wing and wake are 
modeled as a sheet of vorticity. The wing portion, commonly referred 
to as the bound-vortex sheet, is specified, and as a result, a finite 
pressure jump exists across it. The wake portion, commonly referred 
to as the free-vortex sheet, is not specified but rather is force-free 
and is formed as part of the solution in such a way that no pressure 
jump exists across the wake. 

4.2 A Description of the Method 
4.2.1 The wing representation 

The vortex-lattice representation of a typical wing is illus- 
trated in Figure 4.1. The wing is modeled by a lattice of discrete 



► N 
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Figure 4.1. The Unsteady Vortex-Lattice Method 










vortex filaments. Each element of area on the wing is surrounded by a 
closed loop of constant-circulation vortex segments. These segments 
join at the nodes, where the loops make sharp (typically 90°) turns. 

A node is the intersection of the segments and is represented by a 
heavy dot in the figure. Control points are located at the center of 
each element; two typical control points are shown in the figure 
(denoted by small crosses). 

The individual discrete vortex segments connecting the nodes are 
typically members of two different loops (the exceptions are the seg- 
ments along the edges). Consequently, the circulation around an 
individual segment is the algebraic sum of the circulations around the 
two loops to which it belongs. In Figure 4.1 the T's represent 
circulations around the individual segments, and the G's represent 
circulations around the closed loops of segments. For example, 
referring to the figure, Tj is Gg - G^. 

4.2.2 Kinematic considerations 

Two coordinate systems are defined: the fixed inertial frame and 

a moving frame fixed to the wing and aligned with the wing root and 
elastic axis. This would be a deforming coordinate system for the 
elastic-wing model. The problem is posed in terms of the moving ref- 
erence frame. As described in the chapters addressing the equations 
of motion, the dynamic model will be posed in terms of a translating-, 
but nonrotating, coordinate system. 

The position of a point may be described as 
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The velocity of this point is given by 


P 


R + v + 0 ) x r 


(4.2) 


where v is the velocity vector relative to the moving frame and 
is the angular velocity of the moving frame. The components of the 
velocity vectors are expressed in terms of the moving frame. The 
components of the angular velocity in the moving frame may be ex- 
pressed in terms of the derivatives of the Euler angles. A complete 
description of the coordinate transformations and Euler angles for 
general motion is presented in Konstadinopoulos 's work. The trans- 
formations required for this research involve only the pitch angle. 
These transformations are described in the chapter which presents the 
development of the governing equation of motions. 


4.2.3 The Biot-Savart Law 

The equation that serves as the foundation of the UVLM is the 
Biot-Savart Law (see Karamcheti (1980)) which gives the velocity 
induced at point P (see Fig. 4.1) by an individual vortex segment. 


<cos 6 1 


COS © 2 ) 


n » r j 
ft x r. 


(4.3) 


The graphical representation of this equation is illustrated in Fig- 
ure 4.1. T represents the circulation associated with a vortex 
filament described by ft and h is the perpendicular distance from 
the filament to the point P. 


+ 3 
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Equation (4.3) is used to establish the aerodynamic influence 
matrix and the local velocity at a point as affected by all bound and 
free vortex filaments. The velocity field generated identically 
satisfies the continuity equation by this equation. The continuity 
equation for an incompressible fluid is 

V • V = 0 (4.4) 

4.2.4 The boundary conditions 

The flow must satisfy the following boundary conditions. The 
disturbance velocity must approach zero far away from the lifting sur- 
face and wake (satisfied identically through the Biot-Savart Law), and 
the relative velocity normal to the lifting surface must vanish on the 
lifting surface, 

(v - V^ g ) • n * 0 on the wing (4.5) 

where V ls is the velocity of the lifting surface and n is the 
normal to the lifting surface. In addition, for an inviscid fluid the 
Kelvin-Helmholtz theorem, which requires that all vorticity be trans- 
ported with the fluid particles when the pressure is continuous, is 
used to obtain the position of the force-free wake. This theorem may 
be stated as 



where V is the circulation around any arbitrary closed loop. 
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The vortex sheet representing the wing is approximated by a lat- 
tice. Thus, the no-penetration condition, Equation (4.5), is enforced 
at only a finite number of points, the control points of the bound 
lattice. These control points are located at the centroids of the 
elements, and the normal vectors are obtained by forming the cross 
product of the diagonals in each element. 

This lattice serves as a computationally expedient imitation of 
the boundary layers on the lifting surface and the free-shear layers 
in the wake. The velocity field induced by the vorticity is computed 
according to the Biot-Savart law; consequently, the velocity field 
satisfies the continuity equation and decays far from the wing and its 
wake. 


4.2.5 The calculation of the circulation 

The circulations of the loops, the in Figure 4.1, are ob- 

tained from the no-penetration condition 


N 


I A 
j-1 


ij "J 


+ -► 

(V - V ) 
Is w i 


-*■ 

n. 


for i ■ 1,2, ...,N 


(4.7) 


where N is the total number of the elements representing the wing, 
and the Ajj represent the normal component of the velocity at the 
control point of the i— (receiving) element generated by a closed 
loop of unit circulation around the j— (sending) element. At the 
i— control point (V^ g - V w )^ is the velocity of the lifting surface 


minus the velocity induced by the wake elements 
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For a rigid wing, the problem is posed in a moving coordinate 
system attached to a reference configuration of the wing and the 
are constant. For an elastic wing the must be computed at each 

time step. 

4.2.6 The unsteady Kutta condition 

The Kutta condition requires ACp to be zero along the wing tips 
and trailing edge. It is satisfied at each time step by convecting 
the vorticity along the sharp edges of the wing into the wake at the 
local particle velocity. Each node is displaced according to 

Ar = v it (4.8) 

where Ar is the displacement of the node, v is the local velocity 
relative to the moving coordinate system, and At is the time step. 
The local particle velocity is computed from 

v - V - V. - 2 x ? (4.9) 

A 

-► 

where V is the absolute velocity calculated from the Biot-Savart 
Law. To maintain the wake in a force-free position, each node in the 
wake is also displaced according to Equations (4.8) and (4.9). 

4.3 The Impulsive Start of the UVLM 

Initially, the wing is represented by a lattice. There is no 
wake. Before motion begins all circulations are zero. The instant 
motion begins the circulations in the bound vortex lattice change. A 
discrete vortex line forms along the wing tips and trailing edge. The 
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existence of this starting vortex is dictated by the requirement of 
spatial conservation of circulation. This is the vorticity that is 
convected into the wake from the sharp edges of the wing at the parti- 
cle velocity. 

The entire aeroelastic model has been nondimensionalized by de- 
fining a characteristic length, Z c , and characteristic velocity, U,,,. 
From these definitions the nondimensional time step, At, is defined 
as Z c /U„. As much as possible, it is desirable to have uniform ele- 
ments in the lattices. To achieve this end, we choose the character- 
istic length to be the chord of the rectangular elements. At each 
time step the free-vortex system is convected at the local particle 
velocity. Equation (4.8). As a result the streamwise lengths of the 
wake elements are approximately the same as the bound elements. One 
might also note that an increase in the number of chordwise elements 
is accompanied by a corresponding decrease in the actual time step. 

4.4 Graphical Representations of the Flowfield 

In Figure 4.2 a flowfield predicted by the general UVLM is illus- 
trated. The wing-tip vortex system is clearly evident. It strongly 
influences the velocity at the control points along the wing tips. In 
Figure 4.3 the flowfield is shown in different views. In Figure 4.4 
the strength of the vorticity fields of both the bound and free vortex 
sheets is shown. 

4.5 The Aerodynamic Loading 

The aerodynamic loading is determined by calculating the pressure 
jump across each individual element in the bound (wing) lattice. The 
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Figure 4.2 


Wireframe of Wing and Wake Lattice 
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Figure 4.3. Vortex-Lattice Representations 
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Figure 4.4. Wing and Wake Vorticity for Wings of Different 
Aspect Ratios. 
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pressure is calculated from Bernoulli's equation for unsteady flows. 
The pressure difference between the upper and lower surface of element 
i is given by 




(V - R - u> x r)] 


(4.10) 


where AV is the discontinuity in the tangential component caused by 
the local vorticity. The average circulation around the discrete 
vortex segments on parallel edges of an element is considered to be 
the circulation around a sheet of vorticity parallel to and between 
the same edges. The vorticity is concentrated into cores to facili- 
tate computation of the velocity field. The average circulation is 
related to the jump in the tangential component of velocity perpen- 
dicular to these edges, across the thickness of the lifting surface* 
The reader is again referred to Konstadinopoulos (1981). 

The pressure distribution is multiplied by the area of the 
element to provide the elemental force perpendicular to the surface. 
These forces are integrated along the chord and span to provide the 
aerodynamic forces and moments. 



CHAPTER V 


THE NUMERICAL SOLUTION SCHEME 

5.1 Overview 

The equations governing the motion of the structure have been de- 
veloped. In addition, the model that determines the motion-dependent 
aerodynamic loads has been described. The interaction between these 
aerodynamic loads and wing motion presents an interesting challenge 
since the loads are necessary to predict the motion, yet the motion is 
necessary to predict the loads. Therefore, an interacting numerical 
integration scheme has been developed, which determines the motion, 
loads, and wing/wake geometry simultaneously. The method is based on 
Hamming’s Predictor-Corrector Method as described in Carnahan (1969). 

This integration technique coupled with the nondimensionalized 
form of all governing equations results in a reliable arrangement of 
the problem. The characteristic length, £ c , is the chordwise length 
of a lattice element. One time step is represented by the time neces- 
sary for a fluid particle to convect approximately one £ c . The 
predictor-corrector method does not subdivide a time step of integra- 
tion as do other techniques such as Runge-Kutta. Hence, the chordwise 
size of the convected elements remains nearly uniform for all 
calculations of the loads. 

5.2 The First-Order Differential Equations 

The development of the two-degree-of-f reedom model for the rigid 
wing and the multiple degree-of-f reedom model for the elastic wing 
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both result in a set of first-order differential equations as devel- 
oped in previous sections. The set of first-order equations which 
describe the two-degree-of-f reedom motion (see Eqs. (2.15) and (2.16)) 
can be written as 


h * x x [ (K a + C 5 - nCC )r cos(a + <J>) 

K - r cos (a + <f>) 

+ K{ 5 2 r sin( a +<)>)- K y Y - C y n - CC n cos(a +<(>)}] 

l « 5 — — {(C n + K Y)r cos( a + $) 

K - r cos (a + <J>) ^ ^ 

- r 2 £ 2 sin(u + <(>) cos(a + <|>) - (C^? + K^a) 

+ C(n C + r cos 2 (a + $)C )} 
v m T n' J 

where 

Y - n 
a = 5 

The equations of motion for the elastic wing are repeated for 
convenience. 


(5.1) 


(5.2) 


(5.3) 

(5.4) 
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(3.34) 
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In general, we have a system of first-order differential equa- 
tions of the form 

It f l (t ’ V y 2 V 

ar- f 2< c - V y 2 V 

(5.5) 


f „ (t - V Y 2 V 


where n Is 4 for the two-degree-of-f reedom problem and n is twice 
the number of modes for the elastic wing problem* 

5.3 The Basic Predictor-Corrector Method 


In both sets of equations the right-hand side represents the de- 
rivatives of the dependent variables* The general predictor-corrector 
scheme requires the values of the dependent variables at the current 
and three previous time steps* In addition, the derivatives at the 
current and two previous time steps are required. With this informa- 
tion the computation of the state variables for the next time step can 
be calculated* The approach to obtaining the information necessary to 
start the procedure will be described later in this chapter. The 
integration of one time step will be described when all the necessary 


information Is known 
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The state variables at the end of the current time step are pre- 
dicted by the following equation! 


j,i+1.0 j *i~3 3 


+ 4 At(2f 


- f, 


+ 2f , 


) 


j,i "j.i-l ' “j,i-2 
j * If 2, • • • ,n 


(5.6) 


where the subscript i is the index for the time step and the sub- 
script 0 signifies the initial prediction of the state variables. 
The time interval, At, is unity in our formulation due to the non- 
dimensionalization. The local truncation error, e^ £, is estimated 
at the end of each time step and this truncation error is used to 
modify the predicted values to y*. 




7 J.i+1.0 + 


112 




(5.7) 


Then y* is used to obtain the derivatives (f*). Next, the 
"corrector" equation is used to obtain the updated state variables. 




8 [9y j,i “ y j,i-2 
+ 3At(f j >i+ i,M-l + 2f 




- f 




)] 


(5.8) 


The corrector equation is used iteratively (M is the index for 
this iteration) until convergence on the state variables is achieved. 
The derivatives, f* i+1 M-1 , are computed from y^ i+1 M _j. Hence, 
the corrector equation uses the current information. In addition, it 
is important to note that the aerodynamic model is used each time the 
state variables are predicted or corrected. 
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The wake is convected from the position generated at the end of 
the last time step into the new position. This wake, new state vari- 
ables, and wing geometry are used to compute the loads. Therefore, 
the aerodynamic loads are current. 

Once convergence of the state variables and loads has been 
achieved the local truncation error is estimated. 


e j,i+l “ 121 (y j,i+l,M “ y j,i+l,M-l } 


(5.9) 


A final modification is made to the state variables. 


'j,i+l 7 j,i+l,M " *j,i+l 


(5.10) 


5.4 Starting the Predictor-Corrector Method 

The predictor-corrector procedure requires the derivatives of the 
dependent variables from the current and three previous time steps. 

It also requires the derivatives of these variables for the current 
and two previous steps. Carnahan suggests a Runge-Kutta scheme to 
determine these values as a starter for the predictor-corrector 
technique. We elect to use a Taylor series expansion to establish the 
starting values because this approach only requires the aerodynamic 
loads at integral time steps. 

For the nill equation 


dt 


f n (t, y 2 » •••» Y n > 


(5.11) 
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where the right-hand side of Equation (5.11) is the right-hand side of 
Equations (5.1) through (5.4), or (3.34). Then, Equation (5.11) is 
expressed in difference forms between the 1 st and 2 nd time step. 



f (t, Y. 


V 


or 


(t, Y 


r 2’ 


V 


At + Y 


Now recall, 


Y(t + At) - Y(t) + At Y'(t) +|j-Y"(t) + h.o.t. 


(5.12) 


(5.13) 


(5.14) 


where the prime indicates differentiation and h.o.t. represents the 
higher-order terms. Therefore, 


Y ■ Y + At Y’ + Y" + h.o.t. 
ri 3 1^2 ^ 


(5.15) 


or, considering Equation (5.11), Equation (5.15) becomes 


f * At 
n 2 

Y + At f + — Z -, — + h.o.t, 
n 2 n 2 2! 


(5.16) 


Approximating from Y n ^ for 2At and considering relation (5.14) 


( 2 At ) Y^ 

Y (t + 2 At) ■ Y + 2At Y* + ^ 

n 3 n l n l 2! 


+ h.o.t. 


(5.17) 
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or 

Y - Y + 2At f + 2At 2 f’ + h.o.t. (5.18) 

n 3 n ^ n j n ^ 

If we multiply Equation (5.16) by 4 and subtract Equation (5.18) we 
obtain 

3Y - 4Y - Y + 4f At - 2f At + h.o.t. (5.19) 

n n n„ n. n 0 n. 


Next we make the approximation f 1 * f f ; then Equation (5.19) can be 

n 2 n l 

rewritten as 



!'n, 


At + h.o.t 


(5.20) 


In a similar fashion an expression for Y^ is developed 





At + h.o.t. 


(5.21) 


Hence, the starting values for the general predictor— corrector 
procedure are available. 

The integration of the equations of motion by the predictor- 
corrector method is the same for both the rigid wing on the elastic 
foundation and the elastic wing on the rigid support. However, the 
complete integration process includes the interaction with the aero- 
dynamic model; hence, differences may exist in the complete integra- 
tion process. The convection of vorticity from the wing is performed 
in the wing-fixed (not inertial) reference frame. Therefore, the 
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basic difference in the integration process is whether or not the ref- 
erence system is moving and the associated effect on the convection of 
the flowfield at the local velocity. Interestingly, the two-degree- 
of-freedom rigid wing presents a more complicated integration process 
than the multiple-degree-of-f reedom elastic wing since the elastic 
wing is assumed to be cantilevered from a fixed support and the rigid 
wing is attached to a rotating support. The complete integration 
procedures will be described. 

5.5 Integration Method for the Rigid-Wing System 
5.5.1 First-order convection theory 

The integration technique for the equations developed for the 
rigid wing is illustrated by the flowchart in Figure 5.1. First, new 
state variables are predicted. These variables represent both static 
and dynamic contributions. Next, the wake is convected to the new 
force-free position from the position generated at the end of the last 
time step and the aerodynamic loads are computed. Then, the state 
variables are corrected. It is important to note that the aerodynamic 
model is used each time the state variables are predicted or cor- 
rected. Convergence of both the state variables and the aerodynamic 
loads is required. If convergence is not achieved the flowfield is 
convected from that generated (and converged) at the previous time 
step, new aerodynamic loads are determined, and the state variables 
are again corrected. Convergence is checked. If convergence is 
achieved then the conditions for the beginning of the next time step 


are known 
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Figure 5.1. Flowchart of the Integration Scheme Using Lower- 
Order Convection Theory, Rigid-Wing Formulation. 
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5.5.2 Second-order convection theory 

A higher-order convection strategy was examined (see Fig. 5.2). 

In this approach, a two stage corrector process was developed. In 
stage one, the motion was corrected and the aerodynamic loads deter- 
mined for a wake in a fixed position. That is, the position of the 
wake was not recomputed each time the state variables were determined. 
Stage one was repeated to convergence of both the loads and motion. 

In stage two, the wake was convected by an iterative process to a con- 
verged position with fixed state variables. Then, the field points 
for the wake were convected at a velocity based upon the average of 
the corrected variables and the variables from the previous time step. 
The new position of the wake is checked for convergence with the pre- 
vious position* A new two-stage cycle begins: this new wake was used 

for stage one (motion and loads determined) and the converged new 
motion and loads were then used again for stage two (wake position). 
Motion, loads, and the wake were required to converge before the 
integration process continued to the next integration time step. The 
computational time was increased five-fold; however, the predicted 
results showed no discernible differences over the lower-order con- 
vergence scheme. 

5.6 Solution of the Equations for the Elastic Wing 
5*6.1 Static solution 

The solution of the equations for the elastic wing Is handled in 
a somewhat different manner. The static solution of the elastic wing 
is required prior to the start of the dynamic solution. The analysis 
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Figure 5.2. Flowchart of the Integration Scheme Using Second- 
Order Convection Theory, Rigid-Wing Formulation. 
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of the static problem begins with the specification of a nominal angle 
of attack. The corresponding distributed steady aerodynamic loads are 
computed. This load and the weight are then used to calculate the 
bending and torsional deflections. At this point the wing has a new 
shape, but the loads still correspond to the old shape. New loads, 
corresponding to the new shape, are obtained next; then the corre- 
sponding new shape is obtained. The procedure is repeated until 
either both the shape and loads converge or aeroelastic divergence 
(where the aerodynamic forces exceed the elastic restoring forces) 
occurs. The procedure is illustrated in Figure 5.3. The elements 
with dots at their centroids are part of the bound lattice (wing); the 
others are part of the free lattice (wake). 

5.6.2 Dynamic solution 

The integration technique which provides the dynamic solution is 
illustrated in Figure 5.4. Similarities do exist with the rigid wing 
formulation. The wake is convected to the new force-free position 
from the position generated at the end of the last time step. Unlike 
the rigid-wing formulation, this convection is only performed once. 

The new state variables are predicted for the first integration pass 
and are corrected for all subsequent passes. These state variables 
only represent the dynamic contribution. They are added to the static 
contribution and are used through the modal expansions to generate the 
new wing shape. The wake, new state variables, and wing geometry are 
used to compute the loads. Therefore, the aerodynamic loads are 
current. It is important to note that the aerodynamic model is used 
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STATIC 

RESULTS 


Figure 5.3. Flowchart for the Determination of the Elastic 
Wing Deformations Due to Static Loads. 
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Figure 5. A. The Integration Process for the Equations Governing 
the Dynamic Solution of the Elastic Wing. 
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each time the state variables are predicted or corrected. The state 
variables are corrected, new wing geometry is computed, the aero- 
dynamic loads are updated, and convergence of both the state variables 
and aerodynamic loads are checked. If convergence is not achieved the 
scheme returns to the corrector stage. Convergence is required of 
motion and loads at the end of each time step prior to advancing. 



CHAPTER VI 


APPLICATION OF THE MODEL TO SPECIFIC EXAMPLES 

6.1 Introduction 

In this chapter, the aeroelastic model is used with several exam- 
ples to demonstrate the technique. These examples include cases for 
the rigid wing on an elastic support and the elastic wing on a fixed 
support. These examples illustrate the ability of the technique to 
predict the aeroelastic behavior and the unsteady aerodynamic loads of 
a wing. 

6.2 An Example of Aeroelastic Behavior for the Rigid Wing 

First, we consider an example that is similar to one discussed 
by Fung. A large-aspect-ratio wing that can only plunge and pitch 
about the elastic axis is placed in a uniform steady flow. A cross- 
sectional view of this wing is shown in Figure 2.1. 

Fung's model for the aerodynamic loads is limited to two- 
dimensional flow. Our model includes the effect of the wing tips and 
the wake; hence, these predicted results account for three-dimensional 
and unsteady characteristics. The essential physical properties for 
this example are given in Table 6.1 (these are values used by Fung). 

No structural damping is present in this example; therefore, 
damping of the motion results strictly from aerodynamic effects. 

Later, we will demonstrate the additional effect that structural 
damping contributes to the wing response. In addition, the elastic 
axis and the mass axis are coincident; hence, the pitch and plunge 
degrees of freedom are not coupled inertially. The coupling which 
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Table 6.1. Properties of the Rigid Wing Example 


Wing area 
Wing chord 

Elastic axis/mass axis offset 
Elastic axis location (% of chord) 

Mass of wing 

Mass moment of inertia 
Support translational spring stiffness 
Support torsional spring stiffness 
Translational natural frequency 
Torsional natural frequency 
Density of air 


60 ft 2 
60 ft 
0 

50% 


269 


lb sec 
ft 


150630 lb sec 2 ft 
208.5 lb/ft 
363020 ft lb/rad 
.880 rad/sec 
1.552 rad/sec 
.002378 lb sec 2 /ft 4 
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occurs between the degrees of freedom is a result of aerodynamic 
effects* 

Two different freestream speeds are considered for the example 
and the density of air at sea level is used. These conditions are 
chosen to bound the critical dynamic pressure which we define as the 
dynamic pressure at which the motion neither grows nor decays. This 
critical dynamic pressure results in wing flutter. 

Typically, for the two-degree-of-f reedom model, a simulation is 
performed by establishing the flowfield for a wing which is fixed at a 
static angle of attack. Next, the wing is released from this fixed 
position and given an initial disturbance. The response of the wing 
is examined; a decay of motion in both degrees of freedom indicates 
stability and a growth in motion indicates instability. 

6.2.1 Subcritical response without structural damping 

In Figure 6.1 we show the response of the wing to a small initial 
disturbance which has the following form: 

Y(0) - 0 Y(0) - 0.01 

and (6.1) 

<*(0) - 0 a(0) - 0.02 

The freestream speed for this case is small and results in a dy- 
namic pressure which is less than the critical dynamic pressure. The 
displacements and velocities associated with both degrees of freedom 
clearly decay with time. The pitch and plunge displacements reflect 
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the static and dynamic contribution; hence, in time the motion will 
decay and will eventually converge to the static (which includes the 
effects of both aerodynamic and mass induced loading) condition. 

In addition, as one would expect, the frequency of oscillation 
for each degree of freedom closely matches the uncoupled frequencies 
of oscillation associated with the respective elastic springs. The 
aerodynamic normal-force and pitching-moment coefficients are also 
shown. This aerodynamic loading decays to the steady aerodynamic 
loads. 

As a matter of interest, in nondimensional form the frequencies 
are dependent upon the freestream speed; therefore, a change in the 
freestream speed will change the nondimensional frequency of oscil- 
lation. However, for subcritical conditions the frequency in physical 
time would remain nearly constant. Also, referring to the nondimen- 
sional form of the governing equations of motion (see Chapter II, 

Eqs. (2.15) and (2.16)) one finds that the stiffness of the support 
spring is reduced as the freestream speed is increased - in fact, the 
square of the speed appears in the denominator. The nondimensional 
structural damping, if present, is inversely proportional to the 
freestream speed. 

6.2.2 Supercritical response to a small initial disturbance 

Next, we show the response of the wing to a dynamic pressure that 
is larger than the critical dynamic pressure. The initial conditions 
are the same (see Eqs. (6.1)) as those used in the previous case. The 
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time histories of both degrees of freedom and their time derivatives 
are shown in Figure 6.2. 

The initial disturbance grows and, as we will show, can be ex- 
pected to reach a limit cycle eventually. A gradual transition occurs 
very early in the motion and is noticed in the response associated 
with the plunge degree of freedom. The motion for the two degrees of 
freedom changes from the natural frequencies in these modes to a com- 
mon frequency, which is nearly the frequency of flutter in the limit 
cycle. This common frequency lies between the two natural frequencies 
and is closer to the natural frequency in pitch than to the one in 
plunge - a characteristic of the flutter phenomenon which has been 
well documented. The time histories of the aerodynamic normal load 
and pitching moment coefficients are also shown. It must be empha- 
sized that the aerodynamic loads and wing motion are dependent upon 
each other and a form of the solution for the loads or motion is not 
assumed. At each time step, the integration process iterates to a 
converged solution, predicting both the flowfield and the motion of 
the wing simultaneously. 


6.2.3 Supercritical response to a large initial disturbance 

In Figure 6.3, we show the response of the wing to a large 


initial disturbance which has the following form: 


Y(0) - 0 Y(0) - 1.0 

and 


( 6 . 2 ) 


a(0) - 0 


o(0) = 0 
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The dynamic pressure is the same as the previous case. However, 
in this case we introduce a larger initial disturbance to only the 
plunge degree of freedom. 

The equations of motion are uncoupled in the absence of aero- 
dynamic loads since r * 0. Yet, we find that the pitch is strongly 
excited almost immediately, which is clearly the result of coupling 
introduced by the aerodynamic loads. The plunging motion rapidly 
decays. In this case, the transition of the motion to a common fre- 
quency is clearly evident. The delay in this transition (when com- 
pared to the case with the small initial disturbance) is created by 
the large initial condition. 

After the initial rapid decay, the plunging motion appears to 
approach a limit cycle, not to decay to zero. The pitching motion 
continues to grow in time. Presumably, independent of the different 
initial disturbances introduced in these cases, the pitching and 
plunging motion of these two cases will eventually match. In Fig- 
ure 6.4, we show the phase planes for the two cases shown in Fig- 
ures 6.2 and 6.3. If we examine the phase planes associated with the 
plunge motion, we observe that in the upper half of Figure 6.4 the 
small motion appears to be growing around an unstable focus towards a 
limit cycle, while in the lower half of the figure the large motion 
appears to be decaying toward the same limit cycle. 

6.2.4 Subcritlcal response with structural damping 

No structural damping in the support is included in these pre- 
vious results, though aerodynamic damping is modelled by the UVLM. 
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In fact, in some frequency domain solutions, the damping coefficient 
is used as an index to determine the onset of flutter. The present 
method can also account for damping in the support. In addition, the 
present equations of motion can be readily extended to include non- 
linear stiffness and damping for both degrees of freedom thereby 
providing a means to study nonlinear structural behavior. 

In Figure 6.5, we show the response to the same initial condi- 
tions and same dynamic pressure as used in the case reflected in 
Figure 6.2. However, we now introduce structural damping into the 

problem, for the case illustrated C * C * 0.01. Now instead of 

y 

growing, the initial disturbances decay and the wing does not flutter. 
As one might expect, the flutter boundary is quite sensitive to struc- 
tural damping. 

6.2.5 Comparison with theoretical results 

The unstable results discussed above were obtained by using a 
wind speed of 125 feet per second, and the stable results were ob- 
tained using a wind speed of 40 feet per second. The density of air 
at sea level was used and assumed constant. The critical speed lies 
between these two. Suggested methods for extracting the exact 
critical velocity will be discussed later. Fung found the critical 
velocity to be 162 feet per second which is higher than our predicted 
value. The differences in the two approaches are (1) Fung considered 
an infinite aspect ratio, and we consider an aspect ratio of ten which 
introduces finite wing effects; (2) Fung considered a zero static 
angle of attack, and we consider three degrees which introduce 
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Figure 6.5. Response of the Rigid Wing at a Critical Dynamic 
Pressure With Structural Damping Present. 
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nonlinearities; and (3) Fung's aerodynamic analysis is based upon 
two-dimensional, linear theory, and our method is based on the 
unsteady vortex-lattice concept. 

6.3 An Example of Aeroelastic Behavior for the Elastic Wing 

We now turn our attention to results predicted for example cases 
of our elastic-wing model. These results describe the responses to 
initial disturbances for dynamic pressures below and above the flutter 
boundary. 

The physical properties of this elastic wing are given in 
Table 6.2. Although the structural model allows for the spanwise 
variation of all properties, we choose to demonstrate the procedure 
using a wing with constant sectional properties. The elastic axis 
is located at the sectional midchord, and the sectional center of 
mass is located aft of the elastic axis. The angle of attack is 
15 degrees, this angle is defined as the angle of attack prior to the 
introduction of loading due to mass, steady aerodynamic loads, and 
unsteady aerodynamic loads. 

The natural modes for the structure are required for the expan- 
sion of the dependent variables, w and a. This expansion is given 
in Chapter III. The Hunter method, described by Gray (1987), is used 
to determine the natural modes for the elastic wing. This method 
solves the two-point boundary value problem by a transfer matrix 
approach to the finite difference equations. Of course, the modes may 
be derived by any available computational or analytical method. 
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Table 6.2. Properties of the Elastic Wing Example 


Wing aspect ratio 

Wing chord 

Bending stiffness 

Torsional stiffness 

Elastic axis/mass axis offset 

Distributed mass 

Distributed moment of inertia 


10 . 

1.0 ft 

1.2 x 10 5 lb-ft 2 

7.0 x 10 4 lb-ft 2 
.1 ft 

.537 lb-sec 2 /ft 2 
.125 lb-sec 2 
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In this example we use the first and second bending modes and the 
first torsion mode to represent the deformations of the elastic axis 
of the wing. The normalized amplitudes of the mode shapes are tabu- 
lated in Table 6.3. The natural frequencies, scaled to the funda- 
mental frequency in bending, are 1, 5.91, and 4.55. 

As discussed in Chapter V, the solution to the equations for the 
elastic wing consist of both the static and dynamic solutions# Ini- 
tially, an undeformed wing is placed in a steady flow. The deforma- 
tions due to the static loads and weight are computed# If conditions 
are such that the wing will aeroelastically diverge, this will be 
determined during the computation of the static deformations# The 
dynamic response will be computed about the statically deformed shape. 

6.3.1 Subcritical response of the elastic wing 

In Figure 6.6, we show the response of the elastic wing when the 
dynamic pressure is below the critical dynamic pressure# The initial 


disturbance has the 

following form: 


qj - 2.5 

q 2 - -25 

q 3 - .05 

o 

• 

o 

It 

rH 

. cr 

q 2 “ 0.0 

q 3 - 0.0 


The generalized coordinates q^ and q 2 are associated with the 
first two bending modes, q^ is associated with the first torsion 
mode. These generalized coordinates are the time dependent coeffi- 
cients in the expansions for the dependent variables. These coeffi- 
cients give the dynamic contribution only. 
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Table 6.3. Natural Modes Used to Represent the Elastic Wing 


Station 
( 7 . wing) 

Bending 
First Mode 
(uitjj - .044*) 

Bending 
Second Mode 
(u>b 2 - .260*) 

Torsion 
First Mode 
(w tl - .198*) 

0. (root) 

0 . 

0 . 

0 . 

.1 

.0174 

-.0978 

.1563 

.2 

.0650 

-.3001 

.3088 

.3 

.1380 

-.5169 

.4538 

.4 

.2316 

-.6687 

.5875 

.5 

.3413 

-.6975 

.7068 

.6 

.4630 

-.5753 

.8088 

.7 

.5926 

-.3064 

.8908 

.8 

.7269 

.0776 

.9509 

.9 

.8632 

.5287 

.9876 

1.0 (tip) 

1.0 

1.0 

1.0 


♦Frequencies in radians/nondimensional time. 
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The motion decays for all three generalized coordinates to the 
statically deformed shape. In addition, the frequency of oscillation 
associated with the history of the generalized coordinates occurs at 
the natural, uncoupled frequency associated with each respective mode. 

6.3.2 Critical response of the elastic wing 

We now increase the dynamic pressure to that which creates the 
aeroelastic instability. This dynamic pressure was estimated by a 
method which will be discussed later. The initial disturbance has the 
following form: 


q^ ■ .25 ( j 2 * 0.0 q^ * 0.0 

q ^ = 0.0 =* 0.0 q ^ = 0.0 


(6.4) 


These initial conditions are smaller than those used for the 
subcritical case. The motion predicted by the simulations for the 
critical case, which use the initial disturbance expressed in Equa- 
tion (6.3), grew excessively fast due to the large initial deforma- 
tions. Hence, smaller disturbances are used and reflect the same 
instability but at a lower initial rate of growth. 

In Figure 6.7 the response of the wing is represented by the time 
histories of the generalized coordinates. These coordinates are de- 
fined as in the previous case. The motion associated with the first 
bending mode (q^) decays initially; however, after this initial period 
the motion appears to neither grow nor decay. Several harmonics (from 
the different modes) are embedded in this motion. The motion asso- 
ciated with the first torsion mode grows in time. The frequency of 
oscillation for this mode is near the fundamental torsional frequency. 
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Figure 6.7. Response of the Elastic Wing at a Critical 
Dynamic Pressure. 
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The history of <\2 reflects the coupling necessary for flutter. 
Initially, the motion decays at the frequency associated with the 
second bending mode. Then, a shift in frequency and other harmonics 
appear. Around time step 600 one sees that the response (for all 
three coefficients) grows in amplitude. Finally, the frequency of the 
second mode coefficient coalesces with the torsional frequency. A 
flutter mode is created. The motion oscillates about a non-zero mean 
close to (but not coincident with) the static deflection. This drift 
is a nonlinear phenomenon (see Chapter 4 of the text by Nayfeh and 
Mook (1985)). 

6.3.3 Supercritical response of the elastic wing 

In the preceding case we considered a dynamic pressure very 
slightly above the critical dynamic pressure. We now increase the 
dynamic pressure to a larger value. The initial conditions are the 
same as those described by Equation (6.4). 

The response of the wing to these conditions is shown in 
Figure 6.8. An interesting feature is contained in the predicted 
response of the torsion coefficient, q 3 » The torsional contribution 
continues to grow. In fact, the oscillation grows to angles that 
would certainly cause the flow to separate from the surface, violating 
the assumptions used in the UVLM. However, more important is the fact 
that the motion is continually growing. Most likely, this motion 
would be catastrophic to the structure. 

As a footnote, if we were to continue to increase the dynamic 
pressure it is possible that we could return to an aeroelastically 
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stable regime* This characteristic Is illustrated in Figure 1.3. As 
the freestream speed is increased, one or more modes may cross into 
the unstable regime. These modes may return to the stable regime (as 
does the mode 3 root) and the structure is again stable. For obvious 
reasons, the lowest critical condition is of interest in aircraft or 
wind-tunnel model design. 

6.3.4 Comparison with empirical results 

Harris et al. (1963) described a method to approximate the flut- 
ter speed based upon empirical data. The method uses a baseline 
flutter speed index adjusted according to particular characteristics 
(e.g. aspect ratio, mass ratio, sweep angle, mass distribution, etc.) 
of the configuration under study. These adjustments are based upon 
parametric analyses of experimental and analytical studies. 

For the elastic wing, the results reflect a constant freestream 
velocity which keeps the reduced (nondimensional) frequency nearly 
constant for all cases. The density of the fluid is modified which, 
in turn, modifies the dynamic pressure. 

Using the approach described by Harris, the density required for 
the system to become unstable is approximately 0.026 slugs per cubic 
foot. In our calculations, a density of 0.0066 slugs per cubic foot 
was used for the stable (subcritical) case, a density of 0.035 slugs 
per cubic foot for the unstable (supercritical) case, and a density 
of 0.030 slugs per cubic foot for the marginally unstable (critical) 
case. There is very close agreement between the empirical and 
numerically simulated condition for marginal stability. 



86 


6.3.5 Unsteady aerodynamic loads 

The aerodynamic loads and the motion of the wing are computed 
simultaneously and interactively. The primary objective of the 
numerical simulation is to determine the response of the wing and, in 
particular, determine whether or not the wing is aeroelastically 
stable. However, as a consequence, the unsteady aerodynamic loads are 
computed and the history of the pressure distribution on the wing 
could be used to describe the physics of aeroelastic response. Fur- 
ther, the technique could be modified by removing the structural model 
and restricting the wing motion to specified manuevers. 

In Figure 6.9 we show the pressure distribution of the wing at 
three different times during the wing motion. The spanwise pressure 
distributions are shown for several chordwise stations. The values 
of the pressure differences are for the control points of the wing; 
the control points along the wing tips, leading edge, and trailing 
edge are not located directly on the edges of the wing. 

In the upper portion the pressure distribution is shown for the 
wing due to static aerodynamic loads and the resulting static deforma- 
tions. The pressure does decay toward the wing tips; however, a 
slight rise in pressure is indicated near these edges as a result of 
wing-tip vortices. The effect of the second bending mode is clearly 
evident. The pressure distribution along the chord (for a constant 
spanwise station) reflects that distribution one would expect from 
classical theory; the peak pressures are at the leading edge and these 
pressures decay to a near-zero value at the trailing edge as dictated 
by the Kutta condition. The pressure distribution is also shown at 

r ■ Si 
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Figure 6.9. The Pressure Distribution of the Elastic Wing 
at Three Different Times. 
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two times during the motion of the wing. Again, the effect of the 
second bending mode is clearly evident. The effect of the strong 
wing-tip vortex systems created by the rapid motion is particularly 
evident in the bottom part of the figure. 

In Figure 6.10 the (nondiraensional) time history of the pressure 
distribution for individual elements is traced. We show the pressure 
difference across four elements: the wing-tip leading-edge, the wing- 

tip trailing-edge, the wing-root leading-edge, and the wing-root 
trailing-edge elements. The wing-tip elements experience large- 
amplitude motions; hence, large fluctuations occur. In contrast, the 
wing-root elements do not experience large motion because of the 
fixed-wing boundary condition; hence, smaller fluctuations occur. The 
pressure differences for the two trailing-edge elements are near zero 
which substantiates the presence of the Kutta condition in our 
aerodynamic model. 

In Figure 6.11 the (nondimensional) time histories of the aero- 
dynamic normal-force and pitching-moment coefficients are shown for 
the wing-root, mid-span, and wing-tip locations. These coefficients 
are computed from the pressure distributions associated with the 
column of lattice elements located at these spanwise stations. As one 
would expect, the coefficients oscillate at frequencies associated 
with the natural modes. The peaks do not occur at the same time at 


all stations 




Figure 6.10. The Predicted Pressure Distribution for 
Four Selected Elements. 
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6.4 The Identification of Flutter 

We have shown the response of the rigid wing and the elastic 
wing to conditions which result in both aeroelastic stability and 
instability. However, the technique which determines the precise dy- 
namic pressure necessary for the instability has not been discussed. 
Although not within the scope of this research, ultimately we are in- 
terested in determining the lowest possible dynamic pressure at which 
initial disturbances will not decay. Once this exact condition can be 
ascertained then the effects of angle of attack, aspect ratio, or 
other nonlinearities may be examined. 

Currently, we increase the freestream conditions (dynamic pres- 
sure) from a known stable position, simulate the response of the wing 
to an initial disturbance, and examine the response of the wing. At 
some point it is clear that the motion is growing and that a flutter 
mode exists. However, we have only bounded the critical dynamic pres- 
sure. It is difficult to determine the precise dynamic pressure at 
which the wing becomes unstable. In some cases, the motion associated 
with a degree (or degrees) of freedom appears to grow while the motion 
for another degree (or degrees) of freedom appears to decay. In 
addition, several simulations are required. A method is needed to 
identify this critical point with a minimum of required simulations. 

The frequency domain solutions (Desmarais et al., 1978) determine 
this critical point by tracking the eigenvalues (and associated aero- 
dynamic damping) of the system for increasing velocities. We briefly 
describe these methods in Chapter I, the reader is referred to Fig- 
ure 1.3. Hassig (1971) demonstrated a similar approach by examining 
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the damping associated with the assumed harmonic motion. Time domain 
solutions (Dowell et al., 1980) typically use logarithmic decrement 
approaches to determine the flutter boundary. However, these methods 
assume a form of the solution. 

A method to determine the exact flutter boundary as predicted by 
the numerical simulation has not been fully developed; however, two 
promising approaches have been initiated. The first approach is 
based on a parametric identification procedure which is described by 
Konstadinopoulos (1984) and Elzebda (1986). The second approach is 
based on an energy method. The approaches will be briefly described. 

6.4.1 Parametric identification procedures 

Adopting the approach described by Konstadinopoulos, one can 
model the nonlinear motion of the wing by assuming the normal-force 
and pitch-moment coefficients are nonlinear functions of a, Y, and 
a. An equation for each aerodynamic coefficient is formed which 
contains the higher-order terms. A quintic nonlinearity is suspected; 
therefore, the equations contain all 55 terms up to and including the 
fifth-order terms. 

The aerodynamic normal-force coefficient is represented by the 
following equation: 

• . .2 2 *2 

C ■ a.Y + a„a + a_a + a.Y + a c a + a,a + ... 

n 1 2 3 4 5 6 

, *2 »2 , • 2*2 . , 

+ a c/ Y act + a cc Ya a + h.o.t. 

54 55 


(6.5) 
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The aerodynamic pitching-moment coefficient is represented by the 
following equation: 

=■ bjY + b 2 <x + b 3 « + b 4 Y 2 + b 5 a 2 + b 6 o 2 + ... 

+ b 54 Y 2 aa 2 + b 55 Ya 2 a 2 + h.o.t. (6.6) 

We then fit the aerodynamic coefficient versus time curves - such 
as those curves shown in Figures 6.1, 6.2, and 6.3 - using a least- 
squares technique to determine the coefficients a^ and b^. By 
taking off one term at a time we are able to identify those terms 
which have a negligible contribution. The remaining terms provide a 
simple expression for the aerodynamic coefficients. 

The coefficients a^ and b^ will be unique for a particular 
dynamic pressure, static angle of attack, spring constant, etc. 
Konstadinopoulos classified the coefficients as either restoring terms 
or damping terms. By examining the value and, in particular, the sign 
change of the damping terms versus the critical parameter (the angle 
of attack), he was able to identify the angle of attack at which wing- 
rock occurs. 

Our critical parameter is the dynamic pressure and we are inter- 
ested in the minimum value at which initial disturbances do not decay. 
This parameter will be identified by extracting from the damping 
coefficient versus dynamic pressure curve the condition at which the 
sign for the damping changes. 
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Unfortunately, our results have not yielded the same success as 
the work of Konstadinopoulos. Several complications are introduced 
when the method is extended to the aeroelastic application. First, 
we are working with two aerodynamic coefficients. Second, the terms 
of the fifth-order equations (see Eqs. (6.5) and (6.6)) which are 
eliminated are a function of the dynamic pressure; hence, the form of 
the expression is not quite the same for all cases. Third, the use of 
the method for the elastic wing will certainly be complicated by the 
spanwise variation and phase relationships of the aerodynamic 
coefficients. 

Yet, the concept of examining the character of the damping of 
the system is consistent with the frequency and time domain methods 
addressed in the literature. Therefore, the method may yet be the key 
in future determinations of the exact critical condition. 

6.4.2 Identification using the system energy 

An alternate method is suggested to determine the critical condi- 
tion necessary for an instability. This method determines the total 
work being done on the system by the loads during the simulation. The 
method is based on the concept that an aeroelastic instability exists 
when more energy is being extracted from the freestream by the wing, 
then pumped back into it by the wing. Therefore, an instability is 
present when the total energy level of the wing grows in time. 

In Chapter III we develop the equations of motion for the elastic 
wing. The kinetic energy of the wing is expressed by Equation (3.3) 
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and the potential energy of the wing is expressed by Equation (3.4). 
These equations are repeated here for convenience. 

L 1 

T * / ■=■ [m(y)({w + x (y ) ( cos a)o} 

0 * a 

+ {x (y)(sin a)a} 2 ) + I (y)a 2 ] dy (3.3) 

ot eg 

V - J i [EI(y)w" 2 + GJ(y)a’ 2 
0 

- m(y)g[w + x^Cy) sin a]] dy (3.4) 

The total energy is the sum of the kinetic energy and the poten- 
tial energy. The state variables which have previously been predicted 
through a simulation and describe the motion of the wing are substi- 
tuted into the energy expressions. 

Expressions that give the total energy for the rigid wing are 
available from the formulation of the equations of motion described in 
Chapter II. The method is illustrated through an example utilizing 
the rigid-wing model. 

Figure 6.12 describes the time history of the wing in motion for 
both a stable and an unstable case. The predicted motion of the wing 
(Y,Y,a,a) is used in the expression for the total energy. The stable 
case is characterized by an oscillating, yet decaying, level of energy 
and will reach a constant energy level associated with the steady- 
state (i.e. static) conditions. The unstable case is characterized by 
an oscillating, yet growing, level of energy. If the dynamic pressure 




TIME 


Figure 6.12. The Computed Energy of the Rigid Wing for a 
Subcritical and Critical Case. 
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is slightly above the critical condition then the motion may be 
governed by system nonlinearities and the total energy level will 
fluctuate around a nonzero mean. 

The advantage of the method is that now, instead of examining the 
motion associated with each degree of freedom or fitting the aerody- 
namic coefficients with high-order polynomials, a single quantity is 
used. The method applies to both the elastic wing and the rigid wing 
models . 

The disadvantage is that a transition from the stable to the 
unstable case is still difficult to ascertain. One possible method is 
to measure the decrement associated with the decaying oscillation and 
extract a damping term from the energy curve. Negative damping 
indicates the instability. 

6.5 Wing and Wake Graphics 

An advantage of the general technique is that the time domain 
solution coupled with the UVLM provides the capability of graphically 
depicting wing and wake motion. Now, the physics of aeroelastic 
response can be graphically described. 

We have previously demonstrated the capability (Chapter IV) of 
generating graphical representations of the flowfield. We now address 
the generation of movies of the wing and wake motion as predicted by 
the general model. Recently, we (Strganac and Mook (1987), and 
Strganac, Mook, and Mitchum (1987)) presented movies of the wing and 
wake motion. Both stable and unstable behaviors of the wing were 


shown 


98 


In Figure 6.13 we show the vorticity for the wing at several time 
steps. The frames are extracted from the movie. The view is from 
above the wing surface. Near the top of each frame is the leading 
edge of the wing. The brighter regions, typically those regions 
near the leading edge, indicate the more intense concentrations of 
vorticity. As the wing twists and bends in response to the initial 
condition, the vorticity on the wing changes. It is important to 
reiterate that, in turn, the motion is also a function of this 
vorticity through the unsteady aerodynamic loads. 

In Figure 6.14 we show the vorticity and shape of the wing and 
wake at several time steps (which correspond to Fig. 6.13). The 
history of the previous motion is embedded in the wake. The character 
of the vorticity is also indicative of the current motion of the 
wing. The motion of the wing can also be seen in the sequence. 

For purposes of illustration, the physical properties of the 
system are such that for every three cycles of (first mode) bending 
motion there are two cycles of (first mode) torsion motion. The 
geometry of the wing and wake continues to change. It is interesting 
to note that the wing-tip vortex system breaks down and reforms with 
the bending/torsion motion. The aerodynamic loads which are used to 
calculate the motion reflect this behavior. 
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Figure 6.13. The Strength of Vorticity for the Wing at 
Several Time Steps. 
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Figure 6.14. The Predicted Shape and Strength of Vorticity 
for the Wing and Wake at Several Time Steps. 


CHAPTER VII 


CONCLUDING REMARKS 

A numerical model has been described that predicts the steady and 
unsteady aeroelastic behavior of a wing. The method can predict 
steady-state static responses and transient responses to initial dis- 
turbances. The solution of the governing equations is obtained in the 
time domain. A physical interpretation of the response is aided by 
animation of the wing and flowfield histories. 

An aerodynamic model is introduced into the formulation which 
addresses nonlinearities created by static deformations, angles of 
attack, vorticity dominated flows, and unsteady flowfields. The wake 
is computed as part of the solution; hence, the history of the motion 
is included in the technique. Steady and unsteady aerodynamic loads 
are also computed as part of the solution. 

The numerical model has been developed such that the aerodynamic 
and structural models may be independently modified. As a result, 
nonlinear structural models or other aerodynamic models may be 
introduced in the future. Two structural models are demonstrated: 
the rigid wing mounted to a linear elastic support and the elastic 
wing cantilevered from a rigid support. Both formulations treat the 
nonlinear effects of static angle of attack and displacements. 

The equations governing the motion of the structure are coupled 
with the equations governing the motion of the fluid, and the struc- 
ture and fluid are treated as a single dynamical system. The integra- 
tion of the governing equations using the predictor-corrector tech- 
nique requires convergence of the loads and motion at each time step 
of integration. 
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